home *** CD-ROM | disk | FTP | other *** search
/ Aminet 1 (Walnut Creek) / Aminet - June 1993 [Walnut Creek].iso / aminet / misc / sci / ephem_src_4_28.lha / refract.c < prev    next >
C/C++ Source or Header  |  1992-04-17  |  2KB  |  65 lines
  1. #include <stdio.h>
  2. #include <math.h>
  3. #include "astro.h"
  4.  
  5. /* correct the true altitude, ta, for refraction to the apparent altitude, aa,
  6.  * each in radians, given the local atmospheric pressure, pr, in mbars, and
  7.  * the temperature, tr, in degrees C.
  8.  */
  9. refract (pr, tr, ta, aa)
  10. double pr, tr;
  11. double ta;
  12. double *aa;
  13. {
  14.     double r;    /* refraction correction*/
  15.  
  16.         if (ta >= degrad(15.)) {
  17.         /* model for altitudes at least 15 degrees above horizon */
  18.             r = 7.888888e-5*pr/((273+tr)*tan(ta));
  19.     } else if (ta > degrad(-5.)) {
  20.         /* hairier model for altitudes at least -5 and below 15 degrees */
  21.         double a, b, tadeg = raddeg(ta);
  22.         a = ((2e-5*tadeg+1.96e-2)*tadeg+1.594e-1)*pr;
  23.         b = (273+tr)*((8.45e-2*tadeg+5.05e-1)*tadeg+1);
  24.         r = degrad(a/b);
  25.     } else {
  26.         /* do nothing if more than 5 degrees below horizon.
  27.          */
  28.         r = 0;
  29.     }
  30.  
  31.     *aa  =  ta + r;
  32. }
  33.  
  34. /* correct the apparent altitude, aa, for refraction to the true altitude, ta,
  35.  * each in radians, given the local atmospheric pressure, pr, in mbars, and
  36.  * the temperature, tr, in degrees C.
  37.  */
  38. unrefract (pr, tr, aa, ta)
  39. double pr, tr;
  40. double aa;
  41. double *ta;
  42. {
  43.     double err;
  44.     double appar;
  45.     double true;
  46.  
  47.     /* iterative solution: search for the true that refracts to the
  48.      *   given apparent.
  49.      * since refract() is discontinuous at -5 degrees, there is a range
  50.      *   of apparent altitudes between about -4.5 and -5 degrees that are
  51.      *   not invertable (the graph of ap vs. true has a vertical step at
  52.      *   true = -5). thus, the iteration just oscillates if it gets into
  53.      *   this region. if this happens the iteration is forced to abort.
  54.      *   of course, this makes unrefract() discontinuous too.
  55.      */
  56.     true = aa;
  57.     do {
  58.         refract (pr, tr, true, &appar);
  59.         err = appar - aa;
  60.         true -= err;
  61.     } while (fabs(err) >= 1e-6 && true > degrad(-5));
  62.  
  63.     *ta = true;
  64. }
  65.